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Abstract. The question of coarse-graining is ubiquitous in molecular dynamics. In 
this article, we are interested in deriving effective properties for the dynamics of a 
coarse-grained variable £(x), where x describes the configuration of the system in a 
high-dimensional space R™, and £ is a smooth function with value in R (typically a 
reaction coordinate). It is well known that, given a Boltzmann-Gibbs distribution on 
x € R", the equilibrium properties on are completely determined by the free 
energy. On the other hand, the question of the effective dynamics on £(x) is much 
more difficult to address. Starting from an overdamped Langevin equation onxg R™, 
we propose an effective dynamics for £(x) € R using conditional expectations. Using 
entropy methods, we give sufficient conditions for the time marginals of the effective 
dynamics to be close to the original ones. We check numerically on some toy examples 
that these sufficient conditions yield an effective dynamics which accurately reproduces 
the residence times in the potential energy wells. We also discuss the accuracy of the 
effective dynamics in a pathwise sense, and the relevance of the free energy to build a 
coarse-grained dynamics. 
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1. Motivation 



In molecular dynamics, two types of quantities are typically of interest: averages with 
respect to the canonical ensemble (thermodynamic quantities, such as stress, or heat 
capacity), and averages of functionals over paths (dynamic quantities, like viscosity, 
diffusion coefficients or rate constants). In both cases, the question of coarse-graining 
is relevant, in the sense that the considered functionals typically depend only on a few 
variables of the system (collective variables, or reaction coordinates) so that it would be 
interesting to obtain coarse-grained models on these variables. 



1.1. Coarse-graining of thermodynamic quantities 

Computing canonical averages is a standard task in molecular dynamics. For a molecular 
system whose atom positions are described by a vector x G R™, these quantities read 



$(x)<i/i (1) 

where $ : R n — > R is the observable of interest and fi is the Boltzmann-Gibbs measure, 

dfj, = Z' 1 exp(-/3V(x)) dx, (2) 

where V is the potential energy of the system, {3 is proportional to the inverse of the 

system temperature, and Z = / exp(— (3V{x)) dx is a normalizing constant. Typically, 

x represents the position of N three-dimensional particles, hence x G W 1 with n = 3N. 
All the results we prove are also satisfied if x G TP 1 , where T = R/Z denotes the one- 
dimensional torus. 

As mentioned above, observables of interest are often function of only part of the 
variable x. For example, x denotes the positions of all the atoms of a protein and of 
the solvent molecules around, and the quantity of interest is only a particular angle 
between some atoms in the protein, because this angle characterizes the conformation 
of the protein (and thus the potential energy well in which the system is is completely 
determined by the knowledge of this quantity of interest). We thus introduce the so- 
called reaction coordinate 

£ : R n -> R, 

which contains all the information we are interested in 0. Throughout this article, we 
assume that 

Mil ^ is a smooth scalar function such that, 
^ ^ for all x G R n , < m < |Vf(a;)| < M < oo. 

We have supposed that £ is a scalar function. It is not clear to us whether the results 
of this article can be generalized to the case of a multi-dimensional reaction coordinate. 

| In this article, we do not address the difficult question of how to find a good reaction coordinate. 
See for instance [24] for some discussion on that point. 
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To this function £ is naturally associated an effective energy A, called the free 
energy, such that 

d(£ */i) = exp(-pA(z)) dz, (3) 

where £ * // denotes the image of the measure /i by £. In other words, for any test 
function $ : R -»• R, 

/ $(£(x))Z- 1 exp(-pV(x))dx= f $(z) exp(-pA(z))dz. (4) 

JR n Jr 

Expressions of A and its derivative are given below (see Section [27TT) . 

The interpretation of (jlj) is that, when X is distributed according to the Boltzmann 
measure (jSJ), then is distributed according to the measure exp(— /L4(z)) cte. Hence, 
the free energy A is a relevant quantity for computing thermodynamic quantities, namely 
canonical averages. 

In conclusion, the question of coarse-graining thermodynamic quantities amounts 
to computing the free energy, and there are several efficient methods to perform such 
calculations (see for example [6]). There are also interesting questions related to 
computing approximations of the free energy, especially when the number of reaction 
coordinates is large, for example in polymer science, but this is not the subject of this 
article. 



1.2. Coarse-graining of dynamical quantities 

The objective of this work is to address some issues related to the dynamics of the 
system, and how to coarse-grain it. In short, we aim at designing a dynamics that 
approximates the path 1 1— ► £ (X t ) , where £ is the above reaction coordinate. 

To make this question precise, we first have to choose the full dynamics, which will 
be the reference one. In the following, we consider the overdamped Langevin dynamics 
on state space R n (we will discuss this choice below), 

dX t = -VV(X t ) dt + y^P^dWt, X t=0 = X , (5) 

where W t is a standard n-dimensional Brownian motion. Under suitable assumptions 
on V, this dynamics is ergodic with respect to the Boltzmann-Gibbs measure ((21). Hence, 
for /i-almost all initial conditions Xq, 

lim - / $(X t )dt= [ §{x)dii (6) 
T-+00 1 j j Rn 

almost surely. In practice, this convergence is very slow, due to some metastabilities in 
the dynamics: X t samples a given well of the potential energy for a long time, before 
hoping to some other well of V. 

An important dynamical quantity we will consider below is the average residence 
time, that is the mean time that the system spends in a given well, before hoping to 
another one, when it follows the dynamics (JSJ). Typically, the wells are fully described 
through £ (x is in a given well if and only if £(x) is in a given interval), so that these times 



Effective dynamics using conditional expectations 



4 



can be obtained from the knowledge of the time evolution £(X t ), which is expensive to 
compute since it means simulating the full system. 

In this article, our aim is twofold. First, we would like to propose a one- dimensional 
dynamics of the form 

dy t = b{y t ) dt + v / 2^ ZT a(^) dB t , (7) 

where B t is a standard one-dimensional Brownian motion and b and a are scalar 
functions, such that {vt) Q<t<T is a good approximation (in a sense to be made precise 
below) of (£(-Xt)) <t<r- Hence, the dynamics (I7j) can be thought of as a coarse-grained, 
or effective, dynamics for the quantity of interest. A natural requirement is that (J7J) 
preserves equilibrium quantities, i.e. it is ergodic with respect to exp(— (3A(z)) dz, the 
equilibrium measure of £(A), but we typically ask for more than that. For example, we 
would like to be able to recover residence times in the wells from ([7]), hence bypassing 
the expensive simulation of £ (X t ) (see Section @] for some numerical results on that 
quantity). 

Second, we would like to investigate the relation between (JTj) and the coarse-grained 
dynamics 

dy t = -A'(y t ) dt + \[W l dB u (8) 

which is indeed a one-dimensional dynamics, driven by the free energy, and ergodic for 
exp(— (3A(z)) dz. In other words, what is the dynamical content of the free energy? 
This second question stems from the fact that practitioners often look at the free energy 
profile (i. e. the function z i— > A(z)) to get an idea of the dynamics of transition (typically 
the transition time) between one region indexed by the reaction coordinate (say for 
example {x G M n ; < zq}) and another one (for example {x G M n ; > zq}). If 
£(X t ) follows a dynamics which is close to jSJ), then the Transition State Theory says 
that residence times are a function of the free energy barriers [221 [14"]. and then it makes 
sense to look at the free energy to compute some dynamical properties. It is thus often 
assumed that there is some dynamical information in the free energy A. 

The difficulty of the question we address stems from the fact that, in general, 
t — > £(X t ) is not a Markov process: this is a closure problem. A first possibility is 
to try and approximate £(Xt) by a process which has some memory in time, typically 
a generalized Langevin equation (see for instance [H [19], and also [15J ) . A standard 
framework is then the Mori-Zwanzig projection formalism, which is described in details 
in [11] . Note also that, since we are interested in reproducing only some output function 
of X t (namely £(X 4 )), tools from the control theory may be used. Such an idea has been 
followed in II3II6]. 

If a time-scale separation is present in the system, then memory effects may 
be neglected. In the sequel, we make such time-scale separation assumptions (see 
assumptions [H2] and [H3] of Proposition 13.11) . which allow us to approximate £(Xt) 
by a Markov process of the type ([7]). We use the framework of logarithmic Sobolev 
inequalities to write these assumptions. It has the advantage that we do not assume 



Effective dynamics using conditional expectations 



5 



to a priori know how to split x between fast and slow modes, or to split the potential 
energy V between fast and slow terms (otherwise stated, the time scale separation is 
encoded in the constants entering the logarithmic Sobolev inequalities, and not inserted 
a priori in the model). In addition, within this framework, we can handle reaction 
coordinates that are nonlinear functions of x, the natural cartesian coordinates of the 
system (see the numerical simulations reported in Section H]). 

Another possibility is to start from a dynamics which includes an explicit small 
parameter, representing a time scale separation. One may then apply an averaging 
principle (see [15] and the references therein for more details along this idea; see also [26] 
for a comprehensive review of the averaging principle, when applied to deterministic and 
stochastic differential equations). In Section T3.2[ we consider such a case of potential 
energy being the sum of two terms of different stiffness, as an example of application of 
our general result (see the potential energy (1471)). Note that, even if we explicitly insert 
a small parameter in V, our model differs from the one considered in [31], where a small 
parameter appears in the potential energy and in the diffusion coefficient. 

Other strategies are to try and identify fast and slow modes in the dynamics (see 
e.g. [33, 20J), or to postulate a parametric form for the effective dynamics and to identify 
its coefficients by numerical simulation on the complete system [271 36J. 

We finish this section by a discussion of the choice of the full dynamics. We chose 
the overdamped Langevin dynamics (jSJ). Other choices can be made, in particular the 
Langevin dynamics, which is closer to a Hamiltonian dynamics and can also be seen as 
a method to sample the canonical measure (see [5] for a review of sampling methods 
of the canonical ensemble, along with a theoretical and numerical comparison of their 
performances for molecular dynamics). From the analysis standpoint, the dynamics we 
chose is much simpler, since the diffusion is non- degenerate (in contrast to the Langevin 
dynamics, which is an hypoelliptic equation). We do not know whether the theoretical 
results presented in this article (such as Proposition 13.11) can be generalized to the case 
of the Langevin dynamics. From a practical viewpoint, it may be possible to use the 
same strategy starting from the Langevin dynamics to write another low-dimensional 
dynamics. We have not pursued in this direction. As an alternative to continuous 
time processes, one can also model the dynamics of a molecular system by a discrete 
time Markov chain, for instance in a discrete state space, where each state represents a 
different metastable configuration of the system J3TJ[32]. In that setting, the question of 
estimating the accuracy of a coarse-grained dynamics has been addressed in [30], where 
similar bounds as those derived in this article are obtained. 

1.3. Statement of the main results and outline 

We propose a way to derive an effective dynamics of the form (J7J). This defines a process 
(y*)t>o, which we compare with (£{X t ))t>o, where X t satisfies (JSJ). Three quantities can 
be typically considered to estimate the distance between y t and £(X t ) (on the time 
interval [0, T]): 
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• [Dl] pathwise convergence: E(sup te / 0T ) ~Vt\ 2 ), 

• [D2] convergence of the laws of paths: ||£(£(X 4 ) < 4 <t) - £((y t )o<t<T)\\TV, 

• [D3] convergence of time marginals: sup ig ( 0T ) ||£(£(X t )) — C(y t )\\TV- 

In the above estimators, we have arbitrarily chosen to measure distances between 
probability measures by the total variation (TV) norm, but other choices could be 
made. Recall that the total variation of a signed measure v is defined by ||^||tv = 
supj eL oo ||j|| LOC<1 J fdu. If v is a measure on W 1 which has a density with respect to the 
Lebesgue measure, then its total variation is just the L l norm of its density. 

It is clear that a bound in the sense of [Dl] implies a bound in the sense of [D2], 
which implies a bound in the sense of [D3]. Conversely, by the Skorohod theorem, a 
bound in the sense of [D2] implies a bound in the sense of [Dl], for some well chosen 
realizations of W t and B t (the brownian motions in (jSJ) and ((7j)), but this theorem is 
not constructive. The most relevant criterion in practice is [D2]. Indeed, the criterion 
[D3] does not account for the correlations in time of the process, which are important 
to understand its dynamical properties. On the other hand, the pathwise convergence 
criterion [Dl] is too strong: practionners in molecular dynamics are rarely interested in 
the trajectory per se. Moreover, [D2] implies the convergence of the law of escape times 
(hence of residence times in the wells), at least if the escape time is (almost surely) a 
continuous function of paths, which holds under some regularity assumptions (see [3j 
Exercise 3.9.10]). 

Our first objective is to propose, in a general case, some sufficient conditions on 
the reaction coordinate for a bound of type [D3] to be satisfied. We are actually able 
to derive an estimate of the difference between the time marginals which is uniform in 
time. Next, on a toy-model, we investigate, both theoretically and numerically, if these 
conditions are sufficient and necessary for [Dl] and [D2] to hold. 

The article is organized as follows. In Section [2J after introducing some notation 
and recalling some basic relations concerning the free energy, we propose a natural 
coarse-graining procedure, which enables us to obtain an effective dynamics of type ([7]), 
where the functions b and a can easily be computed (see Equations (pM]h (125]) and (|26|) ). 
In Section [3], we prove that the solution y t of the effective dynamics ( 1261) is indeed a good 
approximation of in the sense [D3]. Our argument relies on entropy techniques, 

and is very much inspired by P2J E]- In Section HJ we present some numerical results 
obtained on a simple model, where we compare residence times in the potential energy 
wells as predicted by the reference dynamics (IS) and by the one-dimensional reduced 
dynamics (1261) . Section [5] is dedicated to establishing error estimates in the sense [Dl] 
of pathwise convergence, in a specific case. These estimates are illustrated by numerical 
simulations. 
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2. A "natural" coarse-graining procedure 

2.1. Notation 

We gather here some useful notation and results. Let S 2 be the submanifold of W 1 of 
positions at a fixed value of the reaction coordinate: 

S 2 = {x G M. n ; £(x) = z). 

Let us introduce which is the probability measure \x conditioned at a fixed value of 
the reaction coordinate: 

exp(-/?V)|Ve|" 1 ^ 
dns z = — r , (9) 



/ exp(-/3\/)|V^r 1 rf ( T Sz 



where the measure a^ z is the Lebesgue measure on T, z induced by the Lebesgue measure 
in the ambient Euclidean space R" D E 2 . 

We recall the following expressions for the free energy A and its derivative A', also 
called the mean force (see [7]): 



A(z) = -/TMn (^J Z- 1 exp(-^y) {V^da^j 



(10) 



and 

A'(z) = I Fdvte,, (11) 



where F is the so-called local mean force: 



^S-r^fJU (12) 



|V£| 2 r \\V£ 
In view of f fTOj) . note that (Q reads 

exp(-/j\/)|ver i ^ Sz ,,o, 

Zexp(— pA) 

These expressions can be obtained by the co-area formula [TD] , which we now recall: 
Lemma 2.1 For any smooth function $ : W 1 — > M., 

/ $(x)\V£{x)\dx= / / Qda^dz. (14) 

Remark 2.1 (Co-area formula and conditioning) The co- area formula shows that 
if the random variable X has law ip(x) dx in MJ 1 , then £(X) has law i[>*(z) dz, with 



[ *P\V{\- l da^ 



It also shows that the law of X conditioned to a fixed value z of £(X) is fij] z , where 
is defined by |P|). The measure \VC,\~ 1 da-£ z is sometimes denoted by 8^t x )-z i> n the 
literature. 
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From the co-area formula, we get the following result: 
Lemma 2.2 For any smooth function \ : MJ 1 — > K, consider 



f xiver 1 ^. 



The derivative of reads: 



d X * 



dz 



— + X div 



Proof: For any smooth test function g : 
formula (flUl. that 



L, we obtain, using the co-area 




X|V£| g'(z)daT; z dz = x(x) g'(^{x))dx. 



Hence, 



X , {z)g'{z)dz = I x(x)g'(£(x))dx 
Jm.™ 

= [ xiver 2 ve-v(^oe) 

= - / ^o^div (x|V£|- 2 V£) 

= ~J R 9{z) i div (xi^r a vo ^ 



- / 9W 



s z L 



fa- + X div 



|V£| 



|V£| 



which yields the result. 



□ 



2.2. A non-closed equation 

Consider X t that solves (151). By a simple Ito computation, we have 

d£(X t ) = (-W ■ V£ + /T^) (X t ) dt + |Ve|(X t ) d£? t (15) 

where _B 4 is the one-dimensional Brownian motion 



(16) 



Of course, equation f|T5|) is not closed. Following Gyongy [13] , a simple closing procedure 
is to consider y t solution to 



where 



dy t = b{t, y t ) dt + jW 1 y t ) dB u 

b(t, y)=E [(- W ■ V£ + /T^) (X t ) | = y] 



;i7) 
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and 

a 2 (t,y)=E[\Ve(X t ) | £(**)= v] ■ (19) 

Note that b and a depend on t, since these are expected values conditioned on the fact 
that £(Xt) = y, where the probability distribution function of X t of course depends on 
t. 

As shown in [13] , this procedure is exact from the point of view of time marginals, 
i. e. [D3] in our above classification. This is stated in the following lemma: 

Lemma 2.3 The probability distribution function ip^ of £(X t ), where X t satisfies |3]) ; 
satisfies the Fokker-Planck equation associated to l[T7\): 

d t ipt = d z (-6 + p- l d z {a 2 ^)) . (20) 

Proof: Let us denote ip(t,x) the probability distribution function of X t . It satisfies the 
Fokker-Planck equation 

d t ip = div (VV ip + p^V-ip) . (21) 

In view of Remark 12. 11 the probability distribution function tp^(t, z) of £(Xt) is given by 

it(t,z)= [ ^(^oiver 1 ^,. 

Jt, z 

Using Lemma [2.21 with \ = ')> we obtain 

s^(m) = / e ( vt y +<"v (p^)«V))|Vr 1 ^,.(22) 

By definition, we have the following expressions for b and a in terms of ip: 
Kt,z) =^ r7) J (-VV-Vf+zr'AfJIVfr 1 V<i^„ 

Using again Lemma [2.21 with % = |V£| 2 ^>(£, •), we obtain 

d z {a 2 ^) = d z j |Ve|V^s,= / (Vf- V^ + ^AOIV^- 1 (23) 
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Let us now prove a variational formulation of (I20p . For any test function g, we have 

d f d f 

-77 / ^{t,z)g(z)dz = — / ip(t,x)g(£(x))dx 
at J R at J Rn 

= I div + P^Vip) g o £ dx 

JR n 

= - I (ijVV + /r x V^) • Vf g'oe^x 
jR n 

= - I I \V£\- X tyVV • Ve + Z^W-Ve) da^g'(z)dz 
Jr is, 

= -/T 1 / d z (a 2 ^)g'(z)dz 

+ iV^rM-W- V^ + r'Ae) iPdax z g'{z)dz 

Jr is, 

= -/T 1 / d z (a 2 ^)g'(z)dz + / & g\z)dz. 



This shows that ^ satisfies (12"U1) . □ 
^.5. y4 closed effective dynamics 



The problem with equation ffTTj) is that the functions b and a are very complicated to 
compute, since they involve the full knowledge of ip. Therefore, one cannot consider (jTTl) 
as a reasonable closure. A natural simplification is to consider a time-independent 
approximation of the functions b and a. Considering (ITS]) and (1191) . we introduce (E M 
denoting a mean with respect to the measure /j.) 

b(z) = E, [(- W • V£ + /r x A£) (*) I £P0 = z] , 

= / (-vy-ve+r'Ae)^, (24) 



and 



a 2 (z) 



(|V£| 2 (X) I £(X) = z) , 

|V£| 2 d^ z , (25) 



where /i£ z is defined by ©• This simplification especially makes sense if £,(X t ) is a slow 
variable, that is if the characteristic evolution time of £,(X t ) is much larger than the 
characteristic time needed by X t to sample the manifold S z . This is quantified in the 
sequel. 

In the spirit of ( TlTT) . we next introduce the coarse-grained dynamics 



dyt = b(y t ) dt + o-(y t ) dB t , y t=0 = £(X ). (26) 



The Fokker-Planck equation associated to the above dynamics will be useful. It reads 

d t <p = d z (-b<p + (3' 1 d z (a 2 <P)). (27) 

Let us first prove that the dynamics (l2"6"j) is ergodic for the equilibrium measure 
£ * /i. The distance between y t and £ (Xt) is estimated in Section [31 



da 
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In view of assumption [HI] and of (j25p . we observe that the diffusion coefficient 
of (1261) satisfies a{y) > m > for any y. Hence, the process defined by (126]) is irreducible, 
and admits a unique invariant probability measure. In the following lemma, we prove 
that exp(— (3A(z)) dz is a stationary measure for (1261) . Hence, the process y t defined 
by (1261) is ergodic with respect to this probability (see Has'minskii [18J, Kliemann [21] 
and the references therein). 

Lemma 2.4 The measure £ * fi on M, which has the density exp(— f3A), is a stationary 
measure for If2h]) . 

Proof: We infer from (1251) and ([TBI that 

a 2 exp(-pA) = Z' 1 [ |Vf | exp(-(3V) do^. 

Using Lemma [2.21 with x = Z 1 |V£| 2 exp(— /3V), we obtain 
p-%(a 2 exp(-/L4)) 

= (3~ l Z- 1 ! [V£ ■ V(exp(-(3V)) + exp(-/3V)A£] |Vf I" 1 da Ez , 

</E z 

= z- 1 [ (-vf ■ w + /r x A£) exp(-^y) i ver 1 

is, 

= 6 exp(-/M). (28) 
As a consequence of the above equation, (1271) can be recast as 

d t <t> = d z (-60 + /T% (a 2 exp(-(3A) exp((3A) 0)) 

= /3~% [a 2 <9 z (0exp(/M)) exp(-(3A)] . (29) 
It is now clear that = exp(— (3 A) is a stationary solution of the above equation. □ 
In view of ( l29l) . we observe that = exp(— (3 A) is not only a stationary measure 
for (1261) . but also satisfies a detailed balance condition ((y t ) is a reversible process with 
respect to exp(— f3A{z)) dz). 

Remark 2.2 Let us set f(t,z) = <f)(t, z) exp(f3A(z)) and let g : R — ► R 6e a (time- 
independent) test function. Then a weak formulation of / flffij) zs 

^ / ^)exp(-/3A(0))^ = -/T 1 / a 2 ^ exp(-/M), 

«t </R JR 

which can be rewritten as 

^ / 7(t,e(x))^(a;))exp(-^(a;))^ = -r 1 / V(/o0 • V(#o£) exp(-/3V).(30) 
Jr™ JR" 

The above weak formulation should be compared with the weak formulation of the Fokker- 
Planck equation (21\) associated to $\): 

1 f g exp(-(3V) = -r 1 [ Vf-Vgexp(-(3V), (31) 



dt 

where f = ipexp(/3V), ip is the probability distribution function of X t satisfying (TJp, 
and g : R n — > R is a (time-independent) test function. We observe that is $3l\) for 
functions which depend on x only through £(x). 
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We now discuss the relation between the dynamics (1261) that we propose and the 
dynamics (jSJ). If the function £ is such that |V£| = 1, then a — 1, and in view of (fTTj) . 
(fl2"]) and (Ell), we have 6 = —A'. Hence, in this case, the effective dynamics ( 1261) is 
exactly jSJ). The fact that |V£| = 1 is equivalent to say that £ is the signed distance 
to the submanifold So = {x; £(x) = 0}. Examples of such reaction coordinates include 
C(x u ...,i n ) = x u or f (ar) = \x\. 

More generally, assume that £ is such that a — 1. Then, in view of (1281) . we have 
b = —A', and again ( 1261) is exactly (jBJ). Note however that, in general, a is not a constant 
function, and f|26]) differs from (jSJ). We will confirm in Section H] that (1261) and (jSJ) may 
lead to significantly different numerical results. 

Remark 2.3 Note that a — 1 writes 

[ exp{-/3V) |V£| rfa Sz = / exp(-/3V) | V^l" 1 

Differentiating this equality with respect to z yields ( using again Lemma \2.0j) 
[ (-W ■ V£ + P^AZ) exp(-/?y) | V^- 1 da Sz 

= -jf (^^-rMiv^^ ex P (-w iver 1 *^, 

which is exactly b = —A' . 

Actually, using the fact that £ is a scalar function, it is possible to recover the case 
a = 1 (for which the effective dynamics is driven by the free energy) by two different 
methods. It is not clear to us whether such a reformulation is also possible in the case 
of a multi-dimensional reaction coordinate. 

A first method is to introduce the following reindexation of the foliation (£ z )zeR. 
We set 

h(x) = / o-~\y) dy 
Jo 

and we introduce the new reaction coordinate 

C = ho£. 

Note that the foliation associated with ( is exactly the same as the one associated with £ 
since h : R — > R is a one-to-one function. It is then easy to check that the coarse-grained 
dynamics associated with the reaction coordinate £ is 

dy t = -A'{yt) dt + ,/2fF^dB t , (32) 

where A is the free energy associated to (. We hence obtain a dynamics of the type (jBJ), 
with an appropriate noise (that is, dB t in (j3"2"J) and dW t in (jSJ) are linked by (fTBJ) ). 

Another possibility is to keep £ as the reaction coordinate, and to consider, instead 
of (jSJ), the dynamics 
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The measure // is also invariant for this dynamics. Then, following the same coarse- 
graining procedure, based on the reaction coordinate £, one ends up with the coarse- 
grained dynamics 

dy t = -A'(y t ) dt + jW 1 dBt, 

where A is the free energy associated to £. This is exactly (jSJ), again with an appropriate 
noise. 

3. Error estimation in terms of time marginals 

In this section, we establish conditions on £ under which the effective dynamics ( |26|) is 
close to the dynamics of from the time marginals viewpoint ([D3] in our above 

classification). 

3.1. Error estimation 

Let i/j^(t, z) be the probability distribution function of £(X t ), where X t follows (j3J), and 
4>(t,z) be the probability distribution function of the solution y t to fTSoT) . Our aim is 
to bound the distance, for any time t, between these two one-dimensional probability 
measures. 

We already introduced the total variation norm to measure distances between 
measures. In the case of probability measures, there are two other useful quantities. 
The first one is the relative entropy, which is defined by 

H(u\rj)= [ ln(^)du, 



di] J 

for any two probability measures v and rj such that v is absolutely continuous with 
respect to 77. The relative entropy provides an upper-bound on the total variation norm 
distance, by the Csiszar-Kullback inequality: 



\W-Vhv < V^H(u\ V ). (33) 

The second one is the Wasserstein distance with quadratic cost, which is defined, for 
any two probability measures v and rj with support on a Riemannian manifold E, by 



W(y,T}) = J inf / d s (x,y) 2 dir(x,y). 
y *-en(j/,77) _y ExS 

In the above expression, d%(x,y) denotes the geodesic distance between x and y on E, 



ds(x,y) = inf / 



a(t)\ 2 dt; a E C\[0, 1], E), a(0) = x, a(l) = y 



and Il(z/, 77) denotes the set of coupling probability measures, that is probability measures 
7r on E x E such that their marginals are v and 77: for any test function $, 

$(x) dn(x, y) = / <&(x)dv{x) and / <&{y) dn{x,y) = / $(y) dr)(y). 
SxS is isxE is 



Effective dynamics using conditional expectations 



14 



In the sequel, we will need two functional inequalities, that we now recall pQ: 



Definition 3.1 A probability measure rj satisfies a logarithmic Sobolev inequality with 
a constant p>0if, for any probability measure v , 

H(y\ V ) < Yp 1 ^) 
where the Fisher information I(v\rj) is defined by 

I(u\ V ) = J Vln(j| 



dp. 



Definition 3.2 A probability measure rj satisfies a Talagrand inequality with a constant 
p>0if, for any probability measure v , 



W(u, V )<^-H(u\t,). 

We will also need the following important result (see [25], Theorem 1] and [4]): 

Lemma 3.1 If rj satisfies a logarithmic Sobolev inequality with a constant p > 0, then 
rj satisfies a Talagrand inequality with the same constant p > 0. 

Logarithmic Sobolev inequalities are very useful to prove properties concerning the 
longtime behaviour of solutions to PDEs (e.g. long time convergence of the solution of a 
Fokker-Planck equation to the stationary measure of the corresponding SDE). We refer 
to [U El [35] for more details on this subject. 

We are now in position to present the main result of this section. 

Proposition 3.1 Assume that £ satisfies [HI], and that the conditioned probability 
measures /i£ z , defined by (Tj|), satisfy a logarithmic Sobolev inequality with a constant p 
uniform in z: for any probability measure v on S z which is absolutely continuous with 
respect to the measure we have 

[H2] H(y\to.) < -^/W/i S J. 
Let us also assume that the coupling is bounded in the following sense: 

[H3] k = \\V^F\\l°° < oo, 

where F is the local mean force defined by (TJj|). 

Finally, let us assume that |V£| is close to a constant on the manifold S 2 in the 
following sense: 

iver-^oei 



[H4] A 



a 2 o£ 



< oo. 
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Assume that, at time t = 0, the distribution of the initial conditions of (EP and l&b]) are 
consistent one with each other: ip^(t = 0, •) = <f)(t = 0, ■). Then we have the following 
estimate: for any time t > 0, 

M 2 / 2o2 2\ 

E{t) < TZ3 + -At- (H(m - HMt, •)!**)) , (34) 



4m 2 \ p 2 

where E{t) is the relative entropy of the probability distribution function ipt of£(X t ), 
where X t follows (EJ) ; with respect to the probability distribution function <fi of the solution 
Vt to 

E(t) = H (tf(t, •) •)) = J In (^^) tf{t, z) dz. 

Let us comment on these three assumptions. Assumption [H2] means that pj] z , which 
is a measure on the manifold E z , is easy to sample from. In view of (1341) . the interesting 
case is when p is large, and then assumption [H2] implies that there is no metastability 
in the manifold T, z . This amounts to assuming that the overdamped dynamics with 
respect to /i£ z (which lives on S z ) is well-mixing. Note finally that, in view of (fP3"|) . 
the relative entropy H{v\p^ z ) and the Fisher information I{v\p^ z ) entering assumption 
[H2] read 

Jz z V 1 exp(-/M(z)) J 



and 

I{v\Py, 



V Sz In 



f iv^r 1 ^ 



dv 



exp 

where / is the density of v with respect to the measure |V£| -1 crE z , i.e. f 
and Vs z denotes the surface gradient: 

V Sz = PV, where P(x) = Id ^ (x) 

I v^r 

is the orthogonal projector on the tangent space to S z at point x G S z . 

We now turn to assumption [H3]. Consider first the case when x = (xi,^) G M 2 , 
and £(x) = X\. Then F = V X1 V and Vs z F = V X2 F = Va. ia , a V. Requesting 
that k is small hence amounts to requesting that V XlX2 V is small, where x\ is the 
reaction coordinate direction whereas X2 is the direction in S 2 . We hence ask for the 
coupling of these two directions to be small. In particular, in the case when V(x) = 
^x T Hx for some symmetric positive matrix H G IR nxn and = £(xi, . . . , x n ) = 
(x\, . . . ,x p ) for some p < n, we have that V-£ Z E = if and only if the covariance 
Cov M ((Xi, . . . , X p ) , (Xp+i, . . . , X n )) = 0, where X G R™ is distributed according 
to dp = Z~ l exp(— (3V(x)) dx. Hence [H3] means that the variables (Xi, . . . ,X p ), 
which represent the reaction coordinate directions, are decoupled from the variables 
(X p+ x, . . . , X n ), which represent the directions of E z . 

In Section [3~2| we will consider an explicit example, and compute an estimation of 
p and k in that case, which will help understanding the assumptions [H2] and [H3] . 
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The assumption [H4] is technical. Observe that, if |V£| is a constant number in 
each manifold then A = 0. 



Before proving Proposition 13.11 let us comment on the estimate (13411 . Note first 
that this estimate is uniform in time. The initial conditions for (126]) and (jSJ) are such 
that (p{t = 0,-) = ip^it = 0,-), which explains that E(t = 0) = 0. In the longtime 
limit, the estimate (134p is not optimal since we know that both and if>^ converge to 
£ -k fi (see Lemma [2.4p . This implies that limt_ +o0 E(t) = 0, a property that we prove in 
Corollary 13.11 below. 



To prove Proposition 13. 1[ we will need the following lemma: 

Lemma 3.2 Let ip : R x M" — > R be the probability distribution func- 
tion of X t that solves (TJP- The probability distribution function of £(X t ) is 

ip^(t,z)= / ip(t, •)|V£|~ 1 d(7s z; and satisfies 
Jt, z 

+ , U - * (35, 



where A is the free energy fTW and E is the local mean force (C 
Proof: Using (|22|) . we compute 

exp(-/L4)d*(V> e exp(/L4)) 
= d z ^ + pA'^, 

/" Ve-V(^exp(/3y)) 

+ jf (div (^L) ^ **.+PA'tf> 

which yields (|35|) . □ 
We are now in position to prove Proposition 13.11 

Proof: We know that <fi satisfies the Fokker-Planck equation (1291) . and that ip^ satisfies 
the equation (1201) . Thus, we have: 



d z {-b^ + P~ l d z (a 2 In 
/T 1 J d z [a 2 d z ( ( j ) exp(f3A)) e xp(-f3A)] ^, 



Effective dynamics using conditional expectations 

= - J (-6 ^ + /T 1 ^ 2 ^)) d z ln 

+ I3' 1 J a 2 c>,(0exp(/M))exp(-/L4) d z (^j . 

Using (|23il . we have: 
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(ve- v^e Xp (^y))ex P (- / sy)) iv^r 1 rf^ 2 

Jt, z 

[ (V£-V(^exp(/3V))exp(-/3V)) IV^ 1 da 22 
+ Pb(t,z) ^{t,z). 



Thus, it holds: 



1 J J (Ve-V(Vexp(/3F))exp(-/3F)) |V£r x das, cUn ^ 
+/?" i y a 2 d z (<t>exp(j3A))exp(-pA) d z (^j , 




-r 1 
-r 1 /- 2 o 



V£- V(^exp(/3F)) 



ivei 2 

V£-V(^exp(/^)) 



exp(-/3V) ) (|V£| 2 - a 2 (z)) IV^I" 1 dofc, «9,ln 



+/?~ i y a 2 ^(0exp(/L4))exp(-/3A) «9 Z . 



We next use ( 1351) to get: 



f - ^ 
-/r 1 




a" 



V£- V(^exp(/3V)) 
. |V£| 2 

(exp(-/3A)) «9 2 (^ exp(/M)) - /3 ( A'(z 



exp(-/3V) ) (|Ve| 2 -a 2 (^)) IVO- 1 da Sz «9 2 ln( ^ 



+/?" i y a 2 d z ((f>exp(J3A)) exp(-/3A) <9 2 



-/r 1 




a 2 



V^- V^exp^F)) 

ivei 2 



exp(-/?V)) (|V0 2 -a 2 (^)) iVei^da^^ln^ 



r/>^ d z In ( ^- 



+/3~ i y a 2 exp(-/M)^K- 




-/r 1 



+ / a 2 



<9 2 (0exp(/M)) - «9 2 (^exp(/M)) ( ^ 



|V0 2 



exp(-/?V) (|V£| 2 -a 2 (z)) ivei^da^^ln 



■0^ d z In 



-r 1 / a 2 



<9 z ln 
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'V£- V(V>exp(/3V)) 



18 



+ 



2e 2 



a 2 ! i 



|v£| 2 



-/T 1 1 



£i + ^2 



a 2 



Let us first consider the second term of fl36|) . We write, using [H3], that 

2 



_ J a ^|Vt|-^, 



where ips z is the measure ip(t,x) dx conditioned to £(x) = z: 



2 

z 7 ) 



(37) 



Since satisfies a logarithmic Sobolev inequality (assumption [H2]), it also satisfies 
a Talagrand inequality (see Lemma [3. II) . hence 

W(d^ z ,d^J 2 < -H{d^ z \dp^) < i/(dV E .|d^J. 
P P 

Gathering the above inequality with (1371) . we obtain 

2 



A'(z) 



< — J(#sJ^E 2 



Using [HI], we thus bound the second term of fl36|) : 



a 2 A'(z) 



< 



T 2- 

IV 

P 2 




Vv In 



■0 



In 



exp(-/3F) 



exp(-/3V) 



(38) 



We now bound the first term of (1361) using a Cauchy-Schwarz inequality, [H4] and [HI]: 



vei 2 

Ve- Vln(^exp(/31/)) 



ex P (-(3v)) (|v£i 2 -^)) iver 1 ^ 

V£| 2 -a 2 (z))^ iVer 1 rfa Ez 
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< 



< A 




|V£| 



(iver 



a (z 



)) 



4> iVei" 1 da^ 



a 



2 ■ 




Vf • Vln(^exp(/3F)) 



< A 2 M 2 



VC- Vln(^exp(/3F)) 



^ IV^" 1 do* a 2 



(39) 



We infer from fl36l) and the bounds fl38l) and fl39l) that 



~d7 ~ 2e7 



A 2 M 2 



Vf- Vln(^exp(/3F)) 



/3 M 2 k 2 
2T 2 ~p~ 2 



+ 



ivei 2 

|V Ez In (^exp(/?V))| 2 ^ 



1- 



a 2 



-(f) 



Note that 



|Vln (if> exp(pV))\' 



V£- Vln(^exp(/3y)) 



Using the lower bound on |V£| given by [HI], we hence obtain 

P M 2 k 2 



dE p- 1 \ 2 M 2 

— < 

dt 2ei in 2 

-/r 1 



|Vln(Vexp(/3F))rV + 



1 



a 2 ^ 



d,ln 



2e 2 p 2 



|Vln (^exp^V))! 2 ^ 



We now optimize on ex and e 2 by choosing them such that S\ + e 2 = 2 and 

P' 1 A 2 M 2 /? M 2 k 2 ml ' 2A 2 p 2 

This yields E\ 



2ei m 2 2e 2 p 2 



X 2 p 2 + m 2 (3 2 K 2 



thus 



dE B^M 2 

< 5— 

dt Am 2 
P~ X M 2 
4m 2 
M 2 



2 t m 2 j3 2 K 2 



A 2 + 



P z 
m 2 (3 2 K 2 



\Vhi(ipexp(pV))\ 2 ^, 



4m 2 



2 , m 2 f3 2 K 2 \ d 



P z 



WW, 



dt 



We next integrate this equation between and t and use the fact that -E'(O) = to 
obtain (El. □ 



We now prove a corollary of Proposition 13.11 which strengthens its long-time limit 
behaviour. 

Corollary 3.1 In addition to the assumptions of Proposition l3A[ assume that 
[H5] The measure £ * /i satisfies a logarithmic Sobolev inequality with a constant r. 
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Consider again the probability distribution function ip^ of £(Xt), where X t 
follows and the probability distribution function <fi of the solution y t to l[2o}) . They 
satisfy: 

Vt > 0, •) - </>(t, 0||rv < min (d(t), 2C 2 exp^/T 1 *)) , (40) 

for some positive constant R, where 



/ M 2 / m 2 ll 2 K 2 \ 

Ci(t) = J— 2 \v + — J^J mm - m(t, -m, m 

C 2 = max (v/2^(0(O,-)|^), ^ (^(0, , (42) 
twi/j <i// = exp(— /3A(z)) cfe. 

As a consequence of this corollary, we see that lim^oo ||V^(*> ') — -)||tv = 0. 
Proof: We infer from the Csiszar-Kullback inequality and from the bound fl34|) that 

||^ - 0||tv < V / 2^W) < (43) 
where Ci(i) is given by fj4"Tl) . We also have 

11^ - 0||tv < - + ||0 - V% v , (44) 

where = exp(— (3A(z)) dz is the equilibrium measure Let us first upper-bound 

E CG (t) = H tylff) = [ In ( ^— - ) 0. Using Q22D, we compute 



dE, 



CG 



(9 <^ 111 



(3- 1 [ d z [a 2 d z (<j>exp(J3A)) exp(-/M)] In 
Jr 

-/T 1 [ [a 2 d g {</>exp(/3A)) exp(-/L4)] d z 
Jr 



exp(-/3A) 
In 



exp(-/L4) 



< -m 2 ^ 1 / (f) 
Jr 



d 7 



In 



exp(-/3A) 



where we have used that a 2 > m 2 , which is a consequence of [HI] and ( J25l) . Since 

fi^ satisfies a logarithmic Sobolev inequality with constant r, we infer from the above 
dE 

bound that — ; — < — 2 r m 2 (3~ 1 Eqq. Using a Gronwall lemma, we obtain 
at 

H(<j>\fj?) = E CG {t) < E CG {t = 0) ex V {-2rm 2 13-H) = # (0(0, -)|//) exp(-2 r m 2 /T 1 1), 
and the Csiszar-Kullback inequality then yields 



/i ? || Ty < y/2H Wrf) < ^2H(<P(0,-)\^) exp(-r m 2 /T 1 t). (45) 
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We now turn to the term I — 



TV 



For any function \ '■ 



I, define 



X |V£| da^ z , and observe that 



\x{x)\ dx 




\X\ 



ivei 



d<Tx dz > 



X 



ivei 



which also reads |M 



TV 



> II r 



I TV 



We apply this inequality with x = ip — /i: 



Since and the conditional measures /xs z satisfy a logarithmic Sobolev inequality 
(see [H5] and [H2]), and under assumption [H3], we obtain that the measure /i also 
satisfies a logarithmic Sobolev inequality with some constant R > (see [23]). Hence, 
by a computation similar to the one on Ecg, we obtain 



H{?l>\n) <H(iP{O r )\n) ex P (-2Rp- l t), 



hence 



(46) 



||^ - V% v < V^P^l exp(-/?r 1 t)- 
Gathering (jUJ), (SSJ) and (jlB]) . we obtain 

Wtf - MWv < C 2 exp(-rm 2 p- 1 t)+C 2 exp(-R P' 1 1), 

where C% is defined by (1121) . The proof of [231 Theorem 1.2] shows that < i? < rm 2 . 
The above bound then yields ||^ — 4>\\tv < 2C 2 exp(— Rf3~ l t), which, gathered 
with g3l), yields (gQJ. □ 



Estimation of the upper-bound constants of (3j\ ) in a particular case 

In this section, we give a very formal argument to estimate the constants p and k 
entering the bound ( 1341) . in a specific case. Potential energies in molecular dynamics 
are often the sum of several terms, with different stiffness. For instance, the potential 
energy of an alkane chain, in the United Atom model [29], reads 

V{X) = V2(di, i+1 ) + J2 V 3&) + Yl V ^ + Kon-bondcdPO, 

i i i 

where g^+i is the distance between atoms i and i + 1, Q% is the bond angle made by 
atoms i — 1, i and i + 1, whereas 0, is the dihedral angle defined by the atoms i + j, 
j = —1, . . . , 2. In general, V2 is a much stiffer potential than V3, which is itself much 
stiffer than V4. 

A simple toy-model for such potential energies is 

V £ (X) = V (X) + -q 2 (X), (47) 
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where Vq and g are two scalar-valued functions that do not depend on the small 
parameter e (see Equation and Figure [1] below for a precise example of type (JIT]) ). 
For simplicity, we assume here that the reaction coordinate £ does not depend on e, and 
that it is constant on the manifolds S 2 (in assumption [H4], A = 0). Since the relative 
entropy is always non-negative, the estimate (]3"4l reads 



M 2 R 2 k 2 

E(t)<—^H(^,-)\p £ ). 

4 p\ 

We also assume that the initial condition of (jSJ) is well adapted to the Boltzmann 
measure p £ , in the sense that H(ip(0, -)\p £ ) is upper-bounded by a constant independent 
of e. Thus the above bound reads 

E(t) < c4 

Pe 

for some constant C independent of e. Our aim is to roughly estimate the coefficients 
p £ and k £ in terms of s. 

Since e is small, the Boltzmann measure ([2]) concentrates on the manifold where 
q = 0, and locally looks like a Gaussian measure of variance e around that manifold. 
The same holds for p^ z , that is assumed to satisfy a logarithmic Sobolev inequality 
(assumption [H2]). Hence, we typically have p £ = O (1/e). 

We now compute the local mean force, defined by (fT2"l) : 



ivei 



(^) 



- g lr7 . l9 + r=-^ p div 



£ " |Ve| 2 |V^| 2 r VlVei 2 , 

Recall that ^ does not depend on e. If Vg ■ 7^ 0, then F is of order O (1/e), and so 
is K e . On the contrary, if Vg • V£ = 0, then F is of order 0(1) with respect to e, and so 
is K e . 

Let us summarize our discussion. In the case when Vg ■ V£ = 0, it turns out that 
p E is of order while k £ is of order 1, and the estimate ( 134|) reads 

E(t) < Ce 2 

for some constant C that does not depend on e. Hence, as e decreases to 0, the effective 
dynamics ( |26l) becomes more accurate, in the sense of [D3]. In the case when Vg- V£ 7^ 0, 
both p £ and of order 1/e, and the estimate fl34l) reads 

£?(*) < C 

for some constant C that does not depend on e. So the effective dynamics ( 1261) is not 
particularly accurate. In the next section, we numerically confirm that the criterion 

V£ • Vg = (48) 

has indeed a significant impact on the accuracy of the effective dynamics. 
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4. Numerical results: residence time estimation 

Our aim here is twofold. First, we want to check the accuracy of ( 1261) in a sense related 
to [D2], on a simple system, and also compare this effective dynamics with the coarse- 
grained dynamics (JSJ) based on the free energy. Second, we wish to assess the relevance 
of the criterion ( 1481) . It seems to be an important condition for estimates in the sense of 
[D3] to be meaningful. Is it also a necessary and sufficient condition in order to obtain 
accurate dynamical properties ? 

In the following numerical tests, we focus on the residence times. We have indeed 
already underlined that the characteristic behaviour of the dynamics ([5]) is to sample a 
given well of the potential energy, then suddenly hopes to another basin, and start over. 
Consequently, an important quantity is the residence time that the system spends in 
the well, before going to another one. In this section, we describe a numerical example 
where we have studied such quantities, which contain dynamical information, and are 
related to the estimator [D2]. 

Consider the two-dimensional potential energy 

V £ (x, y) = (x 2 -lf + - e {x 2 + y- l) 2 (49) 

which is of the form fj47|) . with Vo(x,y) = (x 2 — l) 2 and q(x,y) = x 2 + y — 1. For any 
e > 0, the potential V e has two local minima, at (x,y) = (±1, 0), and one saddle point, 
at (x, y) = (0, 1) (see Figure[T]). There are thus two basins, namely {(x, y) G M 2 ; x < 0} 
and {(x,y) G M 2 ; x > 0}. Since V e is an even function of x, the residence times in each 
well are equal to each other. Our aim is to compare the residence time computed when 
the full description of the system is used (that is, we simulate the dynamics (j3J)) with 
the residence time computed from a coarse-grained description, according to f[2"6T) or (jSJ), 
for two different reaction coordinates. 

In the case at hand, a natural reaction coordinate is £i(x, y) = x, since the value of 
£i already gives the information that the system is in the right or the left well. In that 
case, |V£i| = 1, hence the effective dynamics ( |26l) is the same as the dynamics (jHJ), that 
is the dynamics driven by the free energy A% associated to £i. This free energy reads 

A 1 (z) = (z 2 - l) 2 + C(0) (50) 

for some constant C{(3) ensuring that / exp(— f3A\(z)) dz = 1. 

Jr 

Note that V£i • Vg ^ 0. In view of the discussion of the previous section, we do 
not expect the effective dynamics based on £i to be very accurate. 

Consider now the function ^{x, y) = xexp(— 2y), which satisfies V^2 • Vg = 0. We 
expect the effective dynamics ( |26i) . based on £ 2 , to be accurate, at least in the sense 
of the estimator [D3] (time marginals). Here, we want to check its accuracy in terms 
of residence times (and hence in a way related to estimator [D2]). Note that, for this 
reaction coordinate, |V^| is not a constant function, hence the effective dynamics ( 1261 ) 
differs from the dynamics (jSJ) for A = A 2 , the free energy associated to £2- 
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Figure 1. Plot of the double- well potential (|49f . For clarity of the picture, we set 
e = 1. 




-2 -1.5 -1 -0.5 0.5 1 1.5 2 



x 

Figure 2. Crosses: plot of the trajectory X t = (x t ,yt) solution to (JSJ), for the 
parameters e — 0.01 and (5 = 3. Dashed lines: level sets of £2- 

We work with the parameters e = 0.01 and (3 = 3. On Figure [21 we plot the 
trajectory solution to (jSJ), as well as the level sets of £2- We can see that the trajectory 
remains close to the line {(x,y); q(x,y) = x 2 + y — 1 = 0} (since e is small), and that 
the level sets of £ 2 are parallel to Vg, which implies that V^2 is indeed perpendicular 
to Vg. 

With the choice we made for f3 and e, the system is metastable. On Figure [3], we 
plot function of time, where X t = (x t ,y t ) satisfies (jSJ). We clearly see that x t 

remains close to -1 (that is, the system is in the left well) for a long time before hoping 
to the right well. 
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Figure 3. Time evolution 1 *— > Xt, for X t — (xt,yt) solution to f5|), for the parameters 
e = 0.01 and (3 = 3. We clearly see metastability. 




Figure 4. Plot of the functions b, a and A' 2 , for the reaction coordinate £2- Note that 
b and A' 2 are odd functions, whereas a is an even function. Note the large variations 
of A' 2 in the neighbourhood of z = 0. 



The functions b and a, as well as the derivative of the free energy A 2 (respectively 
defined by ff24l) . fl25l) and ffTTT) ) are plotted on Figure HI in the case of the reaction 
coordinate £2- 

Remark 4.1 For all the numerical tests reported in this article, the complete 
dynamics $5\) has been integrated with the Euler-Maruyama scheme 

X j+1 = Xj - AtVViXj) + v /2Atf3~ 1 Gj, 

where, for any j , Gj is a two-dimensional vector, whose coordinates are independent 
and identically distributed (i.i.d.) random variables, distributed according to a normal 
Gaussian law. 
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For the reaction coordinate £ 1; the effective dynamics is (Q|) ; that we have 
numerically simulated with the same algorithm as above. We have used the analytical 
expression ( f50j) of the free energy A±. 

For the reaction coordinate £ 2; the free energy derivative A' 2 and the functions b and 
a have been computed using the algorithm proposed in [7|/. We have chosen to work in the 
interval £2 G [—200; 200], and computed A' 2 , b and a on a grid of size Az = 0.1 (except 
in the interval [—0.3; 0.3], where we used a finer grid of size Az = 5. 10~ 3 , since the 
variations of A' 2 , b and a are larger in the neighbourhood of 0). Values of the functions 
for z in-between points of that grid have been obtained by linear interpolation (see 
Figure 0). We have again used the Euler-Maruyama scheme to numerically integrate 
the dynamics l[2b}) . 

All dynamics have been integrated with the time step At = 10 -4 . 

For the reaction coordinate £j, i = 1,2, the left and the right wells are defined as 
the sets {(x,y) G IR 2 ; £i(x,y) < — £- h } and {(x, y) G R 2 ; £i(x,y) > £- h }, respectively. 
We have chosen the threshold values £* h > and ^ h > such that wells are more or less 
the same for both reaction coordinates. To compute the residence time, we proceeded 
as follow, for both reaction coordinates £1 and £2: 

(i) we first generated 15 000 configurations {(xj,^) G M 2 }i<i<i5oocb distributed 
according to the measure /x, and such that £(xi,yi) belongs to the right well, that 
is i{xi, yi ) > £ th . 

(ii) we next ran the dynamics ([5]) from the initial condition y*), and monitor the 
first time at which the system reaches a point (x(rj), y(rj)) in the left well: 
7i = inf{*; £(x(t),y(t))<-£ th }. 

(iii) from these (Tj)i<j<i 50 oo, we computed an average residence time and a confidence 
interval. These figures are the reference figures. 

(iv) we next consider the initial conditions {£(2^, yA G IR} 1<?;<15000 for the effective 
dynamics. By construction, these configurations are distributed according to the 
equilibrium measure of £ * /z, that is exp(— /3A(z)) dz, and are in the right well. 

(v) from these initial conditions, we run the dynamics (1261) or (jBJ), until the left well is 
reached (y(t) < — £ th when working with (126]) . y(t) < — £ th when working with (jSJ)), 
and compute, as for the complete description, a residence time and its confidence 
interval. 

The results we found for the residence time are gathered in Table [TJ We see that, 
when we work with £2 (which satisfies the condition V£2 • Vq = 0) and with the effective 
dynamics (1261) . we can reproduce the reference residence time (32.5 ± 0.5) within an 
excellent accuracy. If we still use the reaction coordinate £2, but consider as the coarse- 
grained dynamics the dynamics (jSJ) driven by the free energy A 2 , then we obtain results 
that are inconsistent with the reference results given by the complete description of the 
system. 
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Reac. Coord. 




Ref. residence time 


Reduced dyn. type 


CG residence time 




0.13 


32.5 ± 0.5 


Eff. dyn. ([26]) 


32.7 ± 0.5 




0.13 


32.5 ± 0.5 


Dyn. ([8]) 


6.4 ± 0.3 




0.5 


31.6 ± 0.5 


Dyn. (J26D = (JHD 


24.4 ± 0.4 



Table 1. Residence times obtained from the complete description (third column) 
and from the reduced description (last column), for both reaction coordinates (and 
both dynamics (|26|) and ([8|) when applicable). The threshold values — 0.5 and 
^2 h = 0.13) have been adjusted so that the reference residence times for both reaction 
coordinates (31.6 and 32.5, respectively) are almost equal. 



Note also that the results obtained with choosing d as reaction coordinate, which 
is such that V£i • Vg 7^ 0, are inconsistent with the reference results (in that case, the 
effective dynamics (126]) is the same as (JH])). Actually, the coarse-grained dynamics does 
not depend on e (since the free energy Ai does not depend on e), whereas the complete 
description does depend on e. 

5. Pathwise convergence 

In this section, we prove pathwise convergence results between £(X t ), where X t 
solves (111), and y t which solves ( ]26l) . for some potential energies of the type (1471) . On 
these specific examples, we obtain stronger convergence results than in the previous 
sections (namely, convergence in the sense of [Dl] rather than in the sense of [D3] or 
[D2] , as in Sections [3] and Hj) . 

Consider the dynamics ([51) , with the potential energy V £ defined by (14T|) . It reads 

dXI = -S7V {X £ t ) dt - ^V(g 2 )(X t £ ) dt + ^2/FWt, A% = X . (51) 

Note that the initial condition is supposed to not depend on e. The limit of X^ when 
e — > has been identified in [7] : it is a process X t solution of a SDE that we write below 
(see equation (I54")) ). and that is such that q(X t ) = for any t. 

Assume now that Xf G R 2 : then X t belongs to the one-dimensional manifold 

M = {X e R 2 ; q(X) = 0} . (52) 

Assume also that the reaction coordinate £ is such that its restriction £|_A4 on M. is a 
one-to-one map from A4 to some subset of R (that is, £ parameterizes M). In that case, 
it is easy to build a reduced dynamics from (|5TI) . in the limit e — > 0: one first lets e go 
to zero, writes the dynamics of X t , and then makes a one-to-one change of variable to 
write the dynamics in term of £ (X t ) . Our aim is to write conditions under which the 
so-obtained dynamics corresponds to (126]) . which amounts to say that the diagram (1531) 



Effective dynamics using conditional expectations 



28 



is a commutative diagram. 



2D dynamics (I5T]) on Xf 
i 

ltd computation 
i 

Nonclosed dynamics on £(X t e ) 
i 

Conditional expectations 
i 

Dynamics ( 1261) on y\ £ (Xf ) : 

dyf = b £ (yf) dt + a £ (yf)dB t 



pathwise 

convergence 
e^O 



►e-<-0 



ID limit dynamics (l54l) on X 4 

I 
I 
I 

One-to-one change of variable: 

I 
I 

Dynamics ( |60|) on 2^ 



5.1. Limit of (51\) in a pathwise sense 

We now proceed in details. For any X £ .M, let 

Vg ® Vg 



P(X) = Id 



|Vg| s 



be the projector on 7xM, the tangent space to A4 at X. Let us define 

Vg 



n 



|Vg| 



and k = div n. 



Let us now introduce the process X t solution to the equation 

dX t = -P (X t ) V (V + /TMn |Vg|) (X t ) dt - ^nndt + ^W 1 P (X t ) dW t , (54) 

with the same initial condition X t =o = Xo as (l5"TT) . Let us assume that this initial 
condition satisfies g(Xo) = 0, and let us fix a time interval [0,T]. Then (see [7j), under 
some regularity assumptions on q and Vo, there exists a constant C that does not depend 
on e such that 

|2 



sup E|X 

te[o,T] 



X f \ <Ce. 



(55) 



Note also that q (X t ) = for any time t. 

Assume now that there exists a one-to-one map 

X :XGM 2 ^(£(X),g(X)), (56) 

which implies that the manifold M. defined by (1521) can be parameterized by £. Then, 
equation (l54l) is equivalent to the dynamics 

d (e (x t )) = V£ (X t ) ■ dX t + /r 1 P (X t ) : V 2 £ (X t ) dt. 

After some tedious but not difficult computations, we see that the above dynamics can 
be written 



d (e (X t )) = dt (X t ) dt + d 2 (X t ) dt + jW 1 1 Vei (Xt) dB t 



(57) 
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with again dB t = — — \X t ) ■ dW t , and 

lVg-Vw Vg-VVb 1 u 1 Vg r V 2 gVg 

d2= -piwr +u ^r-p K w\ + p u |v# ' (58) 

where we set 

u = V£ • Vg. (59) 

Since X t satisfies the constraint q (X t ) = 0, the dynamics ( 1571) can be rewritten only in 
terms of £ (X t J =: Zt, in the form 

dz t = di(z t ) dt + d 2 (z t ) dt + ^(z t ) dB t (60) 

where, for any z, 

d a (z) = d 1 (x-\z, 0)) , d 2 (z) = d 2 { X ~\z, 0)) , 7(z) = |V£| (x" 1 ^, 0)) . (61) 

5.#. Effective dynamics associated to ( f5ij) using conditional expectations 

We now follow the strategy that we have outlined in Section [2731 Starting from ( I3T1) . 
we first compute the time derivative of £{Xf) by an Ito computation, and next take 
the conditional expectations of the drift and the diffusion terms. We hence obtain the 
effective dynamics ( 1261) . where b £ and a £ (that depend on e since the Gibbs measure \i e 
depends on e) are defined by (T24l) and (1251) and read 

6 e (a) = E Me [(-W 6 • V£ + /T 1 ^) (X) I £(X) = a] 
= [(-W ■ V£ + P^AZ) (X) | £pf ) = a] 

-^E Me [(gVg-V£)(X) | £(X) = a] 

= d?(a) - \ [(qVq ■ V£) (X) | £(X) = a] (62) 

a 2 (a)=E M£ (|V£| 2 (X) |£(X) = a), 

where d\(a) = E Me [(-W ■ V£ + /3 _1 A£) (X) | £(X) = a] . It is easy to check that, for 
any a, we have 

d\{a) = di(a) + O(e) and a £ (a) = 7(a) + O(e), (63) 
where d\ and 7 are defined by (IBTJ) . 

5.5. Sufficient conditions for the pathwise convergence to the effective dynamics fiUfy) 

Let us establish sufficient conditions under which the equation (1601) is equivalent to the 
effective dynamics (12"6"|) . in the limit e — ► 0. We hence request that, in the limit e — *■ 0, 
the dynamics (|26|) and (1601) have the same drift and diffusion coefficients. 

We first see that this condition is satisfied for the diffusion coefficients, in view 
of fl63|) : the diffusion coefficient a s of (1261) converges to 7, the diffusion coefficient of (160|) . 
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We now turn to the drift terms, which is d\ + d 2 in the case of (I60p . and b £ given 
by fl62|) for the effective dynamics ff26l) . In view of (163]) . these drift terms are equal, in 
the limit e — > 0, if and only if 

- d 2 (a) = limj E, s [(g Vg ■ V£) (X) | £(X) = a] . (64) 



In view of ( 158]) and (16T]) . we have 

d 2 ( a ) = rf 2 ( x -i( a , 0)) = d 2a ( X ~\a, 0)) + /T 1 ^ (^(a, 0)) , (65) 

where c?2 a and o?2b do not depend on /3: 

Vg ■ VV^o 

^2a = ^ |Vg|2 , (66) 

Vg-Vw M Vg T V 2 gVg 

d2b= ~lVqT- K Wq-\ +U |V# ' (67) 
On the other hand, we compute, for any a, 

[(? Vg ■ VO (X) | £(X) = «] = E Me [g(X) u(X) | £(X) = a] 



A direct computation shows that 



/ qu d/i E ^ a . 

J So 



[(g Vg ■ VO (X) | OX) = a] = ± 5(a) + 0(e 3 / 2 ), (68) 
where 5(a) does not depend on j3 and reads 

5 a = ^ a, + -^-f Ma, + w (a, — ^ a, , 
ag j(a,0) oq oq 

where 

«(C,?)=«(^ 1 (e,9)), (69) 
K)(0<?) — K)(x _1 (£> an d J — detjac Hence, (1641) reads 

- d 2a ( X - 1 (a, 0)) - ^d 2h {x'\oc, 0)) = i 5(a). (70) 

We want to enforce this relation for any (3. Since d 2a , d 2 b and 5 do not depend on (3, 
this yields 

d 2a (x' 1 (a,0)) = and - rf 2b (x _1 (a, 0)) = 5(a). (71) 

In view of ( 166]) . a sufficient condition for the first relation to hold is 

VaGR, M(x -1 (a,0)) = 0, (72) 

where, we recall, u = V£ • Vg and x is sucn that x(X) = (£(X), g(X)). In what follows, 
we now assume that £ is such that (1721) holds. The second relation of (!7T|) now reads 

VaeM, X|^( x -i( a;0 )) = ^(a,0). (73) 

We have thus proved the following result: 
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Proposition 5.1 Consider the two-dimensional dynamics (T57]). and its one- 
dimensional limit (5~4\ ), when e — > 0. On the other hand, consider the one- dimensional 
effective dynamics (ESP, obtained using conditional expectations, and pass to the limit 
e — > in the drift and diffusion coefficients. 

Under the conditions (72|) and (75]) (where u, x o,nd u are defined by (EE) , (5b]) 



and (E2|) respectively), these two dynamics are the same. In addition, for any T > 0, 
there exists C > and e > such that, for all e < e , we have 

sup E\aX^-y £ t \ 2 <Cs, (74) 
te[o,T] 

where X\ solves 137]) and y\ solves the effective dynamics l[2b}) . 

Proof: We only have to prove the bound fl74]) . We infer from (155]) and assumption [HI] 
that 

sup E\axn-^(X t )\ 2 <Ce. (75) 
te[o,T] 

The drift and diffusion coefficients of the effective dynamics on y\ are b E and o~ £ . In view 
of dS2J), P5, (PD, dTDJ) and flM}, the former satisfies 

6 e (a) = d\{a) - - e E ME [(<? Vg • VO (X) | £(X) = a] 

= dt(a) + 0(e) - (3- l E{a) + 0(y/i) 
= di(a) + rf 2 (a) + 0( v / i). 
In view of ( 163]) . the latter satisfies 

cr e (a) = 7(a) + 0(e). 

Hence, the difference between, on the one hand, the drift and diffusion coefficients of the 
effective dynamics (b £ and a e ) and, on the other hand, the drift and diffusion coefficients 
of the equation (16D]) on z t = £ (X t ) (namely d\ + c?2 and 7), is of order O (y/e). We infer 
from this estimate that, on any bounded time interval, 

sup E\y £ t - £ (X t )\ 2 < Ce. 
te[o,T] 

Gathering that estimate with (175]) yields ([74"]) . □ 

In Sections 13.21 and H] we outlined the condition V£ ■ Vg = as an important 
condition to get a good analytical estimate in the sense of [D3], and good numerical 
results in terms of residence times. If u = V£ • Vg = 0, then conditions (1721 and (1731 
are satisfied, and we also get pathwise convergence (i.e. accuracy in the sense of [Dl]), 
in the simple two-dimensional setting considered in this section. 

Hence, the same condition V£ • Vg = appears, independently of the estimator 
([D3], [D2] or [Dl]) that we choose to measure the accuracy of the effective dynamics. 
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5-4- A standard test-case 

Consider the two-dimensional potential energy 

V £ (x,y) = V (x,y)+ n2( - X J y2 , x e R, y G E, (76) 

where f2 is bounded away from and V does not depend on e, and the associated 
overdamped Langevin equation, which defines the process Xf = (xf,yf). The limit 
dynamics on x\ when e — > is well-known in that case (see for instance [28J): it reads 



dx t = - (d x V (x t , 0) + ^P^] dt + y f W l dB t 



(77) 



which is the overdamped Langevin equation associated to the potential V c g(x) = 
V (x,0) + /T 1 lnft(x). 

We now wish to recover that result within our approach. The potential energy (1761) 
is of the form (1171) . with q(x,y) = Q(x)y. We wish to choose the reaction coordinate 

y) = x. Observe then that u = V£ ■ Vg = Q'(x) y ^ 0. Hence the simple sufficient 
condition u = (see end of Section 15.31) is not satisfied. However, it is easy to see that 
the less demanding conditions ( 1721) and ( 173]) are satisfied. Hence, in the limit e — ► 0, 
the effective dynamics (|26|) is accurate in the sense of pathwise convergence. 

For e > 0, the effective dynamics reads 

d& = 6 6 (&) c/t + v / 2 _ 5 rT rf J B t , (78) 

I n / (a)fi(a)y 2, 



with 6 e (o;) = — ^fj, e ( a ,-) yd x V (a,y) + 2 J. A straightforward computation 

shows that limbJa) = —d x Vo(a,0) — \ . Inserting this relation in (1781) . we 
e^o p\l{ct) 

recover (!77|) . 

Hence, taking the limit £ — > in the effective dynamics that we propose, we recover 
a well-known result. 



5.5. Numerical results on the example fcfity 

In the numerical case considered in Section HJ we showed that the reaction coordinate 
£,2{x,y) = xexp(—2y) satisfies the relation V^2 ■ Vg = 0. In view of Proposition 15.11 
we hence expect good results when working with £2, in terms of pathwise convergence. 
We have checked this as follows. First, we have simulated a solution of (1511) (with the 
potential V e defined by (H9l) ). for a given realization of the two-dimensional noise, with 
e = 0.01. From this trajectory X t (we omit here for clarity the dependence with respect 
to e), we obtain the time evolution £%(Xt), and we can also construct the one- dimensional 
noise (I16p . This noise is next used in the effective dynamics (1261) . We compare both 
trajectories on Figure [5j we observe an excellent agreement over 10 6 time steps (the 
trajectories plotted on Figure [5] have been computed with a time step At = 10" 4 , hence 
T = 100 = 10 6 At). 

In Section HI we also considered the reaction coordinate £i(x, y) = x, which is 
such that Ui(x, y) = V£i ■ Vg = 2x^0. With this choice of reaction coordinate, 
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Figure 5. Comparison of ^{Xt), where X t solves ([5]), and yt solution of (|26|) with 
the reaction coordinate £2- 




I I I I I I I I I I 

2 4 6 8 10 12 14 16 18 

t 



Figure 6. Comparison of £i(X t ), where X t solves ([5]), and yt solution of (|26|) with 
the reaction coordinate £1. 

X -1 (£?<7) — (£)<?+ 1 — £ 2 )> hence wi(x _1 (£,0)) = 2£, so condition ( 1721) is not satisfied. 
We have numerically performed the same comparison with £1 as the one reported above 
for £2- Results are shown on FigureEJ we observe that the complete dynamics (projected 
on the reaction coordinate) and the effective dynamics disagree, as expected. Note also 
the difference in time ranges between Figures [5] and [6] (the former corresponding to a 
time interval 5 times larger than the latter). 

5.6. Numerical results on a three atom molecule 

We conclude this section by considering a system closer to those considered in molecular 
simulation, although we acknowledge that it is still a toy-example. The system is 
made of three two-dimensional particles at position G M 2 , 1 < i < 3 (hence 
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Figure 7. Comparison of 0(X t ), where X t solves (O, and yt solution of (|26|) with the 
reaction coordinate X i— > 9(X). 

X = (r 1; r 2 , r 3 ) G M 6 ), and submitted to the potential 

V(X) = 1 (||n - r 2 || - 4) 2 + ^ (IN - r 3 || - £ ) 2 + ~^P0 - £ ) 2 

= ^ M*) 2 + 93(I) 2 ) + - 9 )\ 

where 9(X) is the angle between the bonds (ri,r 2 ) and (r 2 ,r 3 ), gipT) = ||ri — r 2 || — 
and (?3pT) = ||r 2 — r 3 || — £ . In the above potential, £ is an equilibrium length whereas 
#o is an equilibrium angle. This potential represents stiff bonds between particles 1 
and 2 on the one hand, and 2 and 3 on the other hand, with a softer term depending 
on the three-body angle 9. To remove rigid body motion invariance, we set r 2 = and 
ri • e y = 0. Then it is easy to see that the angle 9{X) satisfies W • Vgi = V6> • Vg 3 = 0, 
and hence seems to be a good reaction coordinate, in view of the several discussions 
above. 

Numerical experiments confirm this belief: choosing this reaction coordinate, we 
considered the effective dynamics ( |26i) . and compared its solution with the time evolution 
9(X t ), where X t solves (J51) . Results are shown on Figure[7](we worked with the numerical 
parameters e = 10~ 3 , £ = 1, 9 = 1.187 and k e = 208): again, we see a good agreement 
between both trajectories. 
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